Assessing the impact of ambient temperature on the risk of hand, foot, and mouth disease in Guangdong, China: New insight from the disease severity and burden

Background The association between the incidence of hand, foot, and mouth disease (HFMD) and ambient temperature has been well documented. Although the severity of symptoms is an important indicator of disease burden and varies significantly across cases, it usually was ignored in previous studies, potentially leading to biased estimates of the health impact of temperature. Methods We estimated the disability-adjusted life year (DALY) by considering the severity of symptoms for each HFMD case reported during 2010–2012 in Guangdong and used distributed lag-nonlinear models to estimate the association between the daily average temperature and daily DALY of HFMD cases at the city-level. We investigated the potential effect modifiers on the pathway between temperature and DALY and pooled city-specific estimates to a provincial association using a meta-regression. The overall impact of temperature was further evaluated by estimates of DALYs that could be attributed to HFMD. Results The overall cumulative effect of daily mean temperature on the DALY of HFMD showed an inverse-U shape, with the maximum effect estimated to be β = 0.0331 (95%CI: 0.0199–0.0463) DALY at 23.8°C. Overall, a total of 6.432 (95%CI: 3.942–8.885) DALYs (attributable fraction = 2.721%, 95%CI: 1.660–3.759%) could be attributed to temperature exposure. All the demographic subgroups had a similar trend as the main analysis, while the magnitude of the peak of the temperature impact tended to be higher among the males, those aged ≥3yrs or from the Pear-River Delta region. Additionally, the impact of temperature on DALY elevated significantly with the increasing population density, per capita GDP, and per capita green space in parks. Conclusions Temperature exposure was associated with increased burden of HFMD nonlinearly, with certain groups such as boys and those from areas with greater population density being more vulnerable.


Introduction
Hand, foot, and mouth disease (HFMD) is a common infection among children under 5 years old [1]. It is a very contagious condition [1], with large outbreaks reported frequently in multiple Asian countries [2][3][4][5]. China has more than 2 million cases being reported annually, making it a country suffering from the heaviest burden of HFMD worldwide, especially in southern regions such as Guangdong [6]. The severity of symptoms is varying substantially across individuals, with most cases characterized by typical and mild clinical manifestations and recovering in a few days without specific medical care [7,8]. In contrast, some cases may suffer from severe complications which may lead to neurological sequelae and even death [9].
Ambient heat exposure was suggested to be one of the most important environmental factors that were associated with an elevated risk of HFMD in numerous previous studies [10,11]. For example, a systematic review on top of 11 studies previously conducted in Asia-Pacific regions suggested an average 5(2-8)% increase in the incidence of HFMD per 1˚C increase in the temperature [12]. Similar findings were also observed in our previous research [13][14][15][16].
However, both literature and our previous findings were based on the reported number of cases solely, without considering the variation in severity of symptoms across cases. The severity of symptoms is an important indicator of the disease burden, with severer cases that may represent a stronger health impact of heat exposure and result in a higher disease burden. DALY is an established indicator of disease burden and has been widely used for the evaluation of the global and regional burden of diseases since the 1990s [17]. It was reported that the healthy life years (e.g., DALY) lost of a severe HFMD case was 2-4 times higher than that of a mild case [18]. Failing to differentiate cases of different severity or capture the true disease burden may lead to biased estimates of the health impact [19].
To address this important knowledge gap, we designed this study to estimate the daily disease burden (i.e. daily DALY) of HFMD at the city level by considering the severity of symptoms; and assess its potential association with the ambient temperature. We also evaluated how the association was modified by individual-level and city-level factors.

Ethics statement
This study was approved by the Institutional Review Board at the School of Public Health, Sun Yat-sen University (Approval number: 2021-018). All patient information of HFMD surveillance data were collected from China Center for Disease Control and Prevention (China CDC). No confidential information was included and all methods were carried out in accordance with relevant guidelines and regulations. All data employed in this study were deidentified prior to analysis and anonymized confidentiality, thus informed consent was not required.

Study population
This study covered the entire province of Guangdong, one of the largest provinces in South China with a total population exceeding 115 million as of 2019 [20]. We included all the HFMD cases reported to the surveillance system operated by China CDC during 2010-2012. Under the Infectious Diseases Act, the system covered both the clinically diagnosed cases which were defined by the National Clinic Guide as patients with the papular or vesicular rash on hands, feet, mouth, or buttocks, and laboratory-confirmed cases as clinically diagnosed cases with laboratory evidence of enterovirus infections [21]. All the cases were mandatorily reported with information including birth date, onset date, ZIP code for the residential address, and the severity (i.e., mild, severe, and fatal cases).

The outcome data
The primary outcome of this study was the disability-adjusted life year (DALY) which was estimated for each case and then summed up across cases reported for each day in each city. DALY is a measure of overall disease burden, expressed as the number of years lost due to illhealth, disability, or early death. Cases of different severity may be suffering from significantly different burdens of disease. Therefore, compared with the incidence or prevalence, DALY better captures the impact of environmental exposures. We calculated DALY with the following formula as used previously [22].
Where D is the severity weight used to measure the impacts of different severity on health. The severity weights for mild, severe, and fatal cases were set as 0.0285, 0.133, and 1 according to the acute episodes of infectious disease from GBD2013 [23]. r = 0.03 is the discount rate to convert the future benefits (e.g., healthy life years) into present-value terms. Β = 0.04 is the parameter from the age-weighting function regressing the wights on age to capture different social roles for the cases at different ages. C = 0.166 is the correction constant in the ageweighting function. a is the onset age of the case. L is the duration of disability (e.g., illness) or time lost due to premature mortality (i.e., death) which was set as 0.0252, 0.0361, and (74.71-a) years for the mild, severe, and fatal cases according to a prior survey [18]. Of which, 74.71 was the annual average life expectancy at birth from 2010 to 2012, in China [24]. For example, the DALYs for the mild, severe, or fatal cases of one year of age are 0.0001, 0.0008, and 33.6529, respectively. Since the DALY estimate for a fatal case was extremely high while the number of fatal cases was very limited for HFMD, we excluded the fatal cases from modeling to avoid irregular fluctuations and improve the reliability of estimation [25]. To assess the impact of fatal cases on the results, we included the fatal cases in the analysis and found that the overall trend was similar, but with very large fluctuations across cities (see S3 Fig).

Exposure and confounders
We obtained the observational data of weather at the monitoring site level from the China Meteorological Data Service Center (http://data.cma.cn/en), with one site for a city. Among the total 21 cities, Chaozhou and Foshan did not have any monitoring sites. Therefore, we assigned weather observations from the nearest sites (in Jieyang and Guangzhou, respectively) to these two cities. The primary exposure in this study was the daily average temperature (˚C) at the city level.
Consistent with previous studies [26], we controlled for multiple confounders including relative humidity (%), day of the week, the indicator of holidays, and the potential time trend of case numbers. Humidity was another important weather factor that was suggested to be significantly associated with the risk of HFMD [27]. The indicators of weekday and holiday were used to fit the potential temporal pattern in the heat exposure and case reporting. Children may be more likely to expose to ambient heat during weekends and holidays while more likely to visit a doctor during workdays when they have symptoms. The time trend was commonly adjusted in time series analysis to remove the potential seasonality and long-term trends in the case series.

Study design and statistical analyses
We used a time series analysis design by matching the outcome, exposure, and confounders at the daily scale for each city. We used a two-stage framework to investigate the association between DALY and temperature.
In the first stage, we used a distributed lag non-linear model (DLNM) with a Gaussian linkfunction to estimate the cumulative change in DALY per Celsius increase in daily mean temperature (β, DALY/˚C). The model was developed for each city and specified as follows: Where DALY was the daily DALY at the city level, following a Gaussian distribution with μ DALY as the expectation. α was the intercept. β was the regression coefficient of temperature. cb(Temperature) was the cross-basis (cb) function used to fit the association between DALY and temperature, as well as the variation of the association over the lag period. The lag was used to fit the incubation period from the infection to the onset of symptoms which usually was 14 days for HFMD [28]. We fit the association and its variation over the lag period with a natural cubic spline (ns) with 4 and 5 degrees of freedom (df), respectively. DOW was the day of the week (i.e., Monday, Tuesday, . . ., Sunday). Holiday was the indicator of holidays (i.e., 1 = holiday, 0 = non-holiday). ns(Time, 4/year) was a natural cubic spline with 4 df per year (semi-annual seasonality) to fit the potential seasonality and long-time trend in the time series.
In this study, we used the Akaike Information Criterion (AIC) to determine the optimal combination of df for the splines from 2-5 and 3-8 for the association and its variation over the lag period, respectively. We used the generalized variance inflation factor (GVIF) [29] to evaluate the collinearity between the independent variables and confirm that the GVIFs for all factors in the model were below 10 [30]. Using Guangzhou, the capital city in the study area, as an example, we observed a GVIF of 1.129 for temperature, 1.201 for humidity, 1.006 for DOW, 1.063 for holiday, and 1.217 for the time trend, suggesting that collinearity was not an issue for this study.
In the second stage, we pooled the city-level temperature-DALY relationships to an overall association for the entire province using a random-effect multivariate meta-regression fitted through the restricted maximum likelihood (REML) [31]. In this stage, we also investigated how the city-level factors may have modified the association between DALY and temperature exposure. Specifically, we included the following potential effect modifiers respectively to a meta-regression model: population density, per capita GDP (gross domestic product), rate of municipal sewage treatment, garbage disposal rate, per capita green space in parks, passenger transporting rate, longitude, latitude, altitude, temperature, relative humidity, air pressure, rainfall, and sunshine hours. The multivariate extension of the likelihood ratio (LR) test was used to test the significance of meta-predictors and differences between models. The residual heterogeneity was measured and tested by the multivariate extension of I 2 statistics and the Cochran Q test [31].
Furthermore, we estimated the DALY of HFMD that could be attributed to the excess exposure to temperature based on the overall DALY-temperature association developed in the second stage. Specifically, we used a forward perspective method [32] as follows: where AD is the attributed DALY; DALY t is the DALY at time t; l = 0, 1,. . ., L, and L is the maximum lag; β temperature,l is the contributions from the exposure (temperature) occurring at time t to the risk at times t+l.
Results were also stratified by individual characteristics including gender, age (<3yrs vs �3yrs), and regions (the Pearl River Delta region, East, West, North, as displayed in S1 Fig). Data cleaning and analyses were conducted using R software 4.1.1. The statistical tests were considered statistically significant when P < 0.05, on two sides.

Data summary
We observed that there were a total of 0.82 million cases reported in the study area during 2010-2012, with a total DALY of 4,451.161 years (Table 1). Most cases were male (64.24%) and under 3 years old (65.60%). We observed similar trends of DALYs and temperature over time (i.e. dates) in Fig 1, both with a major peak in the summer months of the year. However, we observed another minor peak in September following the major one for daily DALY. The Pearson's correlation coefficient between daily DALY and temperature was 0.556 (P<0.001).

The association between DALY and temperature
We observed an inverse-U shape association between DALY and temperature, both for the city-level estimates and the pooled provincial estimates (Fig 2). Overall, the impact on DALY was significantly increasing with the daily average temperature, with the peak of heat effect identified at 23.8˚C. We estimated that each degree increase in temperature at the peak may result in an increase of 0.0331 (95%CI: 0.0199-0.0463) DALYs lost, relative to the reference temperature (i.e., 2.9˚C). Similarly, we estimated an increase of 0.0012 (95%CI: 0.0006, 0.0018) DALYs lost per degree increase in temperature at 3.5˚C (the 2.5 th temperature), 0.0134 (95% CI: 0.0073, 0.0195) at 10.1˚C (25 th ), 0.0250 (95%CI: 0.0153, 0.0348) at 17.5˚C (50 th ), 0.0250 (95%CI: 0.0153, 0.0348) at 24.8˚C (75 th ), and 0.0154 (95%CI: 0.0098, 0.0210) at 31.4˚C (97.5 th ), as described in Table 2.
Similar associations were observed among subgroups by sex, age, and region (Table 2 and Fig 3). The peaks for males (23.8˚C) and females (23.7˚C) were observed at almost the same temperature, while the magnitude of the peak for males was slightly higher than that for females (β: 0.0211 vs 0.0119). The peak among children aged<3yrs (24˚C) was observed at a slightly higher temperature than those aged �3yrs (23.6˚C), while the magnitude of the peak among those aged �3yrs was greater (β: 0.0204 vs 0.0126). The peaks for different regions were observed at different temperatures (23.9˚C for Pearl River Delta, 24.2˚C for East, 24.9˚C for West, and 22.5˚C for North), and the magnitude of the peak in Pearl River Delta (β = 0.0591) was higher than those for the rest of regions (β = 0.01 for East, β = 0.0163 for West, β = 0.0161 for North).

The modification effect of city-level factors
The potential modification effect of city-level factors was reported in Table 3. Except for the passenger transporting rate (P value of LR test = 0.107), all other factors were observed with statistically significant modification effect in the meta-regression model with all P values of LR We further estimated the DALY-temperature relationship with the city-level factors specified at both low (25 th ) and high (75 th ) levels (S2 Fig), and found a significantly elevated impact of temperature on DALY in the areas with a higher population density, higher per capita GDP, or higher per capita green space in parks.

DALY attributed to the temperature exposure
With a maximum lag of 14 days, the attributable DALYs of HFMD due to temperature were estimated for all cases and subgroups by sex, age, and region (Table 4). Overall, a total of 6.432 (95%CI: 3.942-8.885) DALYs lost (attributable fraction, AF = 2.721%, 95%CI: 1.660-3.759%) could be attributed to temperature. Consistent with trends of association, an inverse-U shape relationship of attributable DALYs was observed. Specifically, we estimated the attributable Regarding the intergroup comparison, similar disparities were observed as those in the association analysis (

Discussions
This study has unique advantages over previous research by accounting for the severity of symptoms and calculating the disease burden of HFMD cases. Although the health impact of ambient heat exposure was well documented, to the best of our knowledge, this is the first PLOS NEGLECTED TROPICAL DISEASES study designed to investigate the heat impact from the perspective of disease burden which more accurately captures the health impact than the traditional outcomes solely based on the reported number of cases.
DALY as an important indicator of disease burden was rarely used for infectious diseases. Among the limited studies focusing on the DALY associated with HFMD, Gan et al. (2015) estimated that the loss of DALYs for severe and mild HFMD cases were 3.47 and 1.76 personyear per 1000 persons, respectively [33]. Another multi-country study evaluated the burden of HFMD in eight high-burden countries (including China) in East and Southeast Asia and estimated that HFMD has been causing 96,900 (95% CI: 40,600-259,000) age-weighted DALYs annually [34]. However, these studies did not further link the burden of HFMD to environmental exposures which may have more public health implications.
We found an inverse-U shape relationship between the DALY of HFMD and the ambient temperature, with a peak of DALY identified at about 24˚C. This finding indicated a monotonous increase in the disease burden of HFMD with temperatures under 24˚C. Although limited evidence on the disease burden of HFMD was available for comparison, our finding was consistent with previous findings based on the incidence rate of HFMD. The association between heat exposure and the risk of HFMD has been well documented. For example, one of our previous studies suggested that each 1˚C increase in temperature was associated with an increase of 1.86% (95% CI: 0.92, 2.81%) in the weekly number of HFMD cases among children under 14 years old [35]. The potential mechanism underlying the adverse impact of heat exposure remains unclear. However, the enhanced activity of HFMD associated virus and increased

PLOS NEGLECTED TROPICAL DISEASES
outdoor activities and exposures were considered the most plausible interpretations [14]. When the temperature continued to increase, we observed a declining trend in the estimated DALY of HFMD, which may be associated with reduced outdoor activities due to elevated heat exposure as well as increased use of air conditioning. Similarly, a study from the non- metropolitan areas of the Beijing-Tianjian-Hebei area also reported that the HFMD incidence no longer increased monotonically with the increasing temperature when the temperature exceeded 25˚C [36]. In addition, we observed similar trends for the DALY estimates attributed to the temperature exposure. According to our estimates, the timing (i.e. the temperature under which the maximum DALY was observed) of the DALY peak was generally consistent across subgroups, with the magnitude of the DALY at the peak being slightly greater among the male and those aged �3yrs and from the Pearl River Delta area. These findings were consistent with the previous findings. For example, Huang et al. reported that with the same temperature exposure, the risk of HFMD among the male tended to be greater than that among the female [37], which potentially resulted from the disparity in outdoor exposure frequency between boys and girls. Similarly, compared with children <3yrs, those �3yrs more likely to attend preschools, and have more outdoor activities and social contacts. Wang et al. observed that the RR estimates of temperature on the same-day HFMD incidence among children aged 3-5yrs were higher than that among those <1yr and 1-2yrs [38]. Regarding the geographical variation in the magnitude of peak DALY, our finding of the higher vulnerability among the children in the Pearl-river Delta region relative to other regions was consistent with our previous findings in the same area. We previously observed that the combined impact of temperature and humidity on the HFMD incidence tended to be higher in the Pearl River Delta than non-Pearl River Delta [21]. A potential possibility is that the Pearl River Delta region is much more prosperous, and children in this area might enjoy more outdoor activities with their family or friends during pleasant days. Furthermore, the Pearl-river Delta region was the most populated place in the study area, which may result in a much greater social contact among the population than in other regions. Similar trends across subgroups were observed for the attributable risk of temperature on HFMD DALY.
The geographical, climate, and socioeconomic factors were suggested to be potential effect modifiers for the DALY-temperature associations across different cities. Of those, the model with per capita GDP achieved the lowest residual heterogeneity. Similarly, Xiao et al. reported that in southern China with higher economic conditions, the temperature that corresponds to the highest epidemiological risk tended to be higher than that in northern China [39]. It is expected as the susceptible population from the areas with a higher level of per capita GDP and passenger transporting rate may have more social contacts. On the other side, the developed economy and established infrastructure in the developed areas are more likely to protect the urbanized population from the adverse impact of weather exposures [36]. Additionally, we found that the geographical and climate factors across cities may also modify the temperature-HFMD DALY association. The finding was consistent with a previous finding in the same study area which indicated that geographical factors such as latitude and longitude were important effect modifiers [40].
Some limitations of this study should be acknowledged. First, individual factors such as housing environment, smoking, drinking, and even educational levels may have significantly influenced the disease burden of HFMD. This study analyzed the data at the city level, which may have obscured some factors through the ecological fallacy effect. Second, because the daily number of cases for each community was insufficient for modeling or stratified analysis, all data were pooled at the city level for each day to achieve sufficient statistical power. Third, estimation of the association may be affected by excluding the fatal cases. However, this bias was more likely to underestimate the health impact of temperature towards the null.

Conclusions
Our study contributes to the limited knowledge on the association between temperature and the disease burden of HFMD epidemics. The moderate level of temperature was attributed to the more DALYs. Furthermore, the population density, per capita GDP, and per capita green space in parks may be the potential driven factors for the variation of DALY-temperature associations across different areas.